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Abstract. - We study drag-induced diffusion of massive particles in scale- free velocity fields, 
where superdiffusive behavior emerges due to the scale-free size distribution of the vortices of the 
underlying velocity field. The results show qualitative resemblance to what is observed in fluid 
systems, namely the diffusive exponent for the mean square separation of pairs of particles and 
the preferential concentration of the particles, both as a function of the response time. 



Introduction. — Understanding how complex flows in their distribution. 



j^Vdvect massive particles is a challenging problem not only 
O when addressing natural phenomena like rain drop forma- 

1 l ion pQ and protoplanetary disks [2], but also for environ- 

^qihental and industrial applications such as dispersion of 
^pollutants [3] and design of the combustion chamber in 

^j- engines [4]. 

During the advection of massive particles by the flow, 
different transport mechanisms come into play such as the 
• dissipative dynamics of trajectory attr actors and the ejec- 
^^tion of the particles from the vortices (eddies) by cen- 
^^trifugal forces. The interplay of these mechanisms, when 
driven by the multi-scale vortical dynamics of turbulent 
^ flows, leads to classical empirical features such as super dif- 
• i-H fusion [5] and preferential concentration (PC) [6]. 

. Experimental and computational studies show that ad- 
C^vection of particles by turbulent flows is generally a su- 
perdiffusive process, that is, the mean squared separation 
of the particles (R 2 ) grows faster than t, the first power 
of time [7]. An early empirical study by Richardson [8] 
suggested that the mean squared separation of elements 
of carrier fluid itself grows with t 3 . 

However, preferential concentrations only occur when 
the density of particles differs from that of the carrier flow. 
The regions of high vorticity act on the heavier particles 
like centrifuges pushing them toward regions of low vor- 
ticity and high strains while trapping lighter particles [6] 
which leads to clustering of particles and inhomogeneity 



Experimental difficulties in measuring Lagrangian 
quantities and limitations of the numerical simulations in 
resolving a wide range of scales of flow motion have posed 
major obstacles in study of these processes [9]. During the 
last decades, great progresses, experimental and computa- 
tional, have been made in overcoming these obstacles (see 
Ref. pUH] for two recent reviews). 

Attempts have been made in using toy models to sepa- 
rate and study the mechanisms present in real turbulent 
flows. Bee et al. [12] introduced a simple flow model to 
study the role of ejection of particles by eddies in the dis- 
tribution of massive particles. The model consists of a grid 
of two-dimensional cells with high vorticities which eject 
particles with a rate depending on the response time of 
the particles to the flow. Their results for the mass prob- 
ability distribution in the fully mixed state show quali- 
tative agreement with what is observed in real systems. 
However, they do not address the problem of existence 
of superdiffusive regimes. Moreover, the model lacks a 
multi-scale vortical structure, a dominant feature in real 
turbulent flows, which indeed turns out to be responsible 
for superdiffusive behaviors. 

In this Letter, we introduce a heuristic model for two 
dimensional scale-free velocity fields which incorporates a 
multi-scale vortical structure, to study the role of ejection 
of particles by eddies, separated from other mechanisms 
such as dissipation. Here, we show that the multi-scale 
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Fig. 1: (Color Online), (a) Modelling multi-scale velocity 
fields, using rigid discs to represent vortices (eddies) of dif- 
ferent sizes. Each eddy has a fixed center and rotates clock- 
wise (green) or counter-clockwise, in order to ensure absence of 
sliding at the contact points, (b) Trajectories of three massive 
particles when particles follow Stokes drag (see Eq. JT])). 

structure is a key feature responsible for superdiffusion 
regimes in the transport of massive particles. In addition, 
within the limitations of the model, the results obtained 
show qualitative resemblance to what is observed in real 
systems. In particular, we discuss the applicability of the 
results to the problem of preferential concentration and 
Richardson's law. 

Bearings as models for velocity fields. — Our 

model for the velocity field of the flow is based on ran- 
dom polydisperse packings of discs fulfilling the particular 
condition of all discs being able to rotate simultaneously 
without sliding on each other at their contact points. At 
their contact points all discs have the same tangential ve- 
locity in order not to have frustration. Therefore, is a 
global parameter of the system while the angular velocity 
of disc i is inversely proportional to its size, = ^iR%- 

Figure [TJl represents such a random packing where col- 
ors indicate the sign of the angular velocity of the discs. 
One notes that no two discs which are in contact will have 
the same color since they rotate in opposite directions. 
Such bichromatic packings are referred to as bearings and 
were originally introduced in the field of granular materi- 
als as a geometrical explanation of seismic gaps in tectonic 
plates [IMS]. 

Random bearings used here are constructed by directly 
imposing the bichromacy condition. Initially a number 
of discs of random sizes are distributed randomly in the 
system. To each disc one color or the other is assigned 
randomly. The rest of the space is filled by inserting the 
largest possible discs with colors chosen such that no two 
discs of the same color touch each other . 

The configurations are complete (space-filling), self- 
similar and contain some level of randomness introduced 
by the initialization and the choice of parameters in the 
construction procedure (see Refs. [T41IT6] for more details.) 
Following Mandelbrot [17 one can define the Hausdorff- 
dimension of the void space, which is related to the ex- 



ponent of the size distribution of the discs. In practice, 
a lower cut-off on the size of the vortices is inevitable, 
leading to a limited range of scales defined as rM/r m . 

In this context, we regard each rotating disc as a vortex 
with a linear velocity profile inside, that is, the fluid in 
each disc is moving in the same way as if it were a rigid 
body. In this way, we obtain a velocity field containing 
eddies of many different scales. Although, the statistical 
properties of such a flow field do not necessarily coincide 
with those of real turbulent flows, its energy spectrum 
follows a power-law whose exponent is directly related to 
the exponent of the size distribution of the eddies (see 
Ref. [13]). During their rotation, the centers of the eddies 
are kept fixed which makes the velocity field stationary. 

The particle-flow interaction depends on the shape of 
the particles and the difference between their density and 
that of the fluid [18]. Here, we will consider particles much 
denser than the carrier fluid. Neglecting buoyancy and 
considering small particle Reynolds numbers, when the 
size of the particles is much smaller than any active scale 
of the flow, the force exerted on the particles by the flow 
follows Stokes drag: 

dV P _ Vf ~ % m 

dt ~ T ( ' 

where r is the response time of the particles and Vf and 
V p are the velocities of the flow and the particle, respec- 
tively. For a single vortex we have Vf = wz x (r — f* c ), 
where r and f c denote the position of the particle and the 
center of the vortex, respectively, and z is the unit vector 
perpendicular to the vortex. Equation (pQ) is solved ana- 
lytically to obtain particle trajectories. Figure [TJd shows 
three typical trajectories of particles through the system 
of eddies. 

Inside the voids between eddies, we assume in Eq. (pp) 
two different cases, namely, ballistic (Vf = V p ) and damp- 
ing (Vf =0). By assuming ballistic trajectories in the 
voids, we exclude all dissipation in the flow field. In the 
damped case the particles are slowed down which corre- 
sponds to a loss of energy and represents viscous dissipa- 
tion on the scale of the smallest eddy, i.e. the Reynolds 
cut-off. 

As the control parameter for flow-particle interaction, 
the Stokes number S is often used instead of the response 
time r. It is dimensionless and related to r through the 
following relation: 

S=-, (2) 

where r c is a characteristic time of the system. However, 
since there is no unique characteristic time in our system 
we use r in our simulations for parameterizing the flow- 
particle interaction. 

Superdiffusion and preferential concentration. — 

Having set up the velocity field as described above, a total 
number of 10 4 particles are released in a small region in 
the center of the system with very small initial velocity in 
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Fig. 2: Different stages during the diffusion of 10 4 particles: 
(a) r = 14, t = 2, (b) r = 14, t = 6, (c) r = 14, t = 10. (d) 
r = 0.01, t = 3. The velocity field is defined through 1.6 x 10 4 
vortices embedded in a circular space of unit size (shown in 
the background) where the centers of the vortex stay fixed in 
time. The size of the vortices varies from tm = 0.02 down to 
r m = 10 -4 in units of system size (rM/r m = 200). 



randomly chosen directions. Figure [2ji, [2Jd, and [2^ repre- 
sent, respectively, the early, intermediate and final stages 
for one specific configuration with a moderate value of r. 
Figure [2]i shows the final stage for a very small value of r 
and here the preferential concentration becomes visible. 

To improve the statistics, we average the results over 50 
different eddy configurations, where the cut-off of the eddy 
size is r rn = 10 -4 , and the maximum size is tm = 0.02 in 
units of the system diameter. In analogy to turbulence, 
the range of spatial scales between the sizes of the small- 
est and the largest discs will be here referred to as inertial 
range. In the simulations, the particles which reach the 
boundary of the system are removed permanently. In or- 
der to prevent this from affecting the results, we gauran- 
tee reaching the mixed state before a signifacant number 
of particles has left the system by choosing the size of the 
largest eddies to be much smaller than the system size. 

Figure [3] shows the mean square pair distance (R 2 ) be- 
tween two particles as function of time. The initial regime 
(dotted line) for very short times, t < to, corresponds to 
the period when the particles have not yet penetrated into 
the eddies of different scales. For large times, t >t m , par- 
ticles either escape or spread uniformly throughout the 
system yielding (R 2 ) ~ L 2 , as indicated by the horizon- 
tal dotted line, where L = 1 is the radius of the system. 
In the intermediate regimes, to < t < t m , where parti- 
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Fig. 3: Mean square separation (R 2 ) between two particles as 
function of time t for a specific response time r = 3. Distances 
are in units of the system size L, the radius of the circular space, 
and time is in units of L/vo, where vo is the tangent velocity of 
the discs. The two times to and t m delimit the interval within 
which particles are in the inertial regime (see text). Time to 
marks the beginning of the inertial regime whereas beyond t m 
particles are uniformly distributed in the system ({R 2 ) ~ 1). 
The dashed line has slope 7 and its intersection with the dotted 
lines enables one to derive to and t m (see text). 



cles are still in the inertial range, one observes (R 2 ) ~ t 7 , 
since these are regimes where the multi-scale structure of 
the flow is dominant. These regimes are limited to times 
when the particles have not yet gone through scales much 
larger than the size of the largest eddies. This behavior is 
qualitatively similar to that of observed for separation of 
fluid particles via direct numerical simulations Ref. [T9] . 

The values for 7, to and t m depend on the response 
time r, and are found using a power-law fit to (R 2 ) in the 
inertial range [to, t m ], which can be written in the following 
form, 



(3) 



where {R 2 )o denotes the mean square separation of the 
particles before entering the inertial range. 

As shown in Figs.H^ andHjD, to and t m behave as power- 
laws of r, namely to = Br^ and t m = Ar a respectively, 
independently of having dissipative voids or vacuum. The 
fit beyond r (dashed lines) yields, A = 4.7 ±1.2, a = 
0.44 ± 0.14, B = 0.437 ± 0.005 and f3 = 0.303 ± 0.005. 
The value for \og(R 2 )o = 7 log (to /tm) ^ 10 -6 does not 
depend on r but only on the specific initialization of the 
system. Thus, from Eq. (j3]) we can write 7 as a function 
of r, 

- a (4) 
1 ~ 6 + (51ogr [ } 

with a = log (R 2 )oj b = log (B/A) and 5 = f3 — a. Figure 
plots 7 as a function of r (symbols) comparing it with 
Eq. (jl]) (dashed line) which fits the observed values of 7 
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Fig. 4: (a) The transient time to as a function of r together 
with (b) the mixing time t m and (c) the exponent 7. Fits in 
(a) and (b) yield power-laws, namely to ~ and t m ~ r a 
(see text). In (c) the dashed line indicates 7 in Eq. @ where 
a = -6, 6 = 2.38 and 5 = 0.101. 
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Fig. 5: (a) The correlation dimension D2 as a function of time 
for three values of r. The inset shows how D2 is determined: 
at each instant t, one computes the average number C(l) of 
pairs of particles that are closer than a distance /; this yields 
a power-law C(l) ~ l D<2 . At intermediate times D2 attains a 
minimum (dashed line) which indicates maximal preferential 
concentration (see text), (b) The correlation dimension 

Dff c) when maximum preferential concentration is attained in 
the inert ial range. The dashed curve fitting D^ c ^ in the range 
r\ < r < T2 resembles the behavior found with DNS calcu- 
lations of particles in homogeneous turbulent flows [20] (see 
text). 



in the middle range of the r-spectrum, approximately at 
T\ < t < T2. The deviations from Eq. (|4]) imply that to 
and tm obey a power-law only for the intermediate and 
not for the limiting values of r which is what one should 
expect; (i) when r — )• 00, particles' inertia is very large 
and therefore they travel ballistically through the system 
without being significantly deflected by the flow, yielding 
7—^2, (ii) when r — > 0, which characterizes tracers (fluid 
particles), the motion of the particles follows Richardson's 
law (R 2 ) oc t 3 0. Although the original Richardson's 
experiment was in three dimensions it can be shown that 
Richardson's law is independent of the dimensionality of 
the system [19]. 

These observations raise the hypothesis that the main 
ingredient underlying Richardson's law is the multi-scale 
structure of vortices with the dissipative dynamics playing 
a minor role. However, this cannot be rigorously verified 
due to the limitations of the model. As mentioned before, 
our velocity fields are stationary and lack the sweeping 
effect. Therefore, the particles become trapped in local 
eddies in the limit of r = and cannot disperse. 

In the second part of our study, we provide an estimate 
of the range of r, (ti,T2), showing that the effect of the 
limitations of the model is minor. This interval is also 
shown in Fig. St. 

Let us next examine the occurrence and properties of 
preferential concentrations. To quantify the degree of pref- 
erential concentration we use the correlation dimension as 



also used by other authors [20 , defined as [2T| [22]: 

d(Io g C(0) 



£>2 = lim 

l->0 



d(logi) 



(5) 



where C(l) is the fraction of pairs of particles separated 
by a distance smaller than I. This quantity takes on the 
value zero, one and two when particles are all clustered in 
a single point, on a line, and uniformly distributed over 
the space, respectively. As can be seen from the inset of 
Fig. Et, in our case C(l) ~ Z^ 2 , which defines an effective 
exponent D2 for each instant t. 

Figure [5^ shows the evolution of D2 in time for three 
values of r. For initial times, the correlation dimension 
depends strongly on the initial condition. For very large 
times, one observes that D2 converges to the dimension of 
the system itself (D2 = 2), i.e. particles distribute homo- 
geneously throughout the system. For intermediate times, 
one observes a local minimum (dashed line), which implies 
that preferential concentration is maximal. For larger sys- 
tems, we observe that the plateau where D2 = be- 
comes broader (not shown). Therefore, we regard = D^f^ 
as the characteristic value of the correlation dimension. 
While the results shown here are for the case of empty (bal- 
listic) voids, similar results are obtained when the voids 
are dissipative. 

Figure [SJd shows as a function of r. Within the 

range t\ < r < T2 one finds a behavior (dashed line) 
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consistent to what is known for the correlation dimension 
in numerical experiments of heavy particles in turbulent 
fluids, namely for r = and r — ^ oo the correlation di- 
mension converges to the dimension of the system (see 
Ref. [20]). 

The deviations from such behavior observed outside the 
range t\ < r < T2 are due to the finite range of scales of 
eddies. As mentioned before, the particles are initialized 
inside a small neighborhood moving in a randomly chosen 
direction. For large r values (r > T2), particle paths are 
approximately ballistic, moving apart from each other as 
a growing circle at large times which yields D2 ~ 1. For 
r < Ti , particle are approximately tracers and get quickly 
trapped in a circular trajectory inside one vortex. This 
last behavior is strengthened when voids are dissipative 
(bullets), since particles are trapped in the voids. 

Conclusions. — We introduced a two-dimensional 
heuristic model for multi-scale fluid velocity fields and 
studied the dependence of the drag-induced diffusion of 
massive particles on it. The model is based on discs of 
different sizes representing flow eddies which roll on each 
other without sliding. By computing the exponent of the 
mean squared displacement as a power-law of time we 
find superdiffusive regimes within a broad range of re- 
sponse times which is consistent with scaling analysis and 
with the particular limiting cases where Richardson's law 
holds and also with the ballistic regime for heavy parti- 
cles. The results show that the main feature responsible 
for the emergence of superdiffusion is the multi-scale vor- 
tical structure of the flow. We also find the characteristic 
footprint of preferential concentration as observed when 
particle disperse in turbulent flows, like for instance the 
plankton in the ocean. Being able to incorporate multi- 
scale features and at the same time separate distinct mech- 
anisms for heavy particle dynamics in fluids, this model 
can be further extend to more general situations where, 
e.g. dissipation occurs inside the vortices. 
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